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Methods and Systems for Tracking of Amplitudes , Phases and 
Frequencies of a Multi- component Sinusoidal Signal 

Related Applications 

This application claims the benefit of U.S. 
5 Provisional Application No. 60/433,808 filed December 17, 
2002 . 

Field of the Invention 

The present invention relates generally to the 
decomposition of an observed signal into its constituent 
10 components and, in particular, to the estimation and 

tracking of characteristic parameters of the constituent 
components of the observed signal . 

Background of the Invention 

The estimation of the characteristic parameters of 
15 a dominant sinusoidal component of a received signal is a 
common feature within radio receivers and many other 
systems. The characteristic parameters of the dominant 
sinusoidal component are typically required to facilitate 
the demodulation and extraction of information carried by 
2 0 the received signal. 

In radio communications Amplitude Modulation (AM) 
and Frequency Modulation (FM) are commonly used, either in 
conjunction with one another or independently. If either FM 
or AM is used exclusively, it is frequency or complex 
25 amplitude information, respectively, that is to be extracted 
from the received signal that has been corrupted after 
traversing a transmission medium (i.e. a channel). Thus, a 
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single independent demodulation method for either FM or AM 
would suffice. 

However, a more challenging problem is to devise a 
method to estimate both the complex amplitude (amplitude and 
5 phase) and the frequency of multiple dominant sinusoidal 
components of a signal, or to devise a radio receiver that 
can, for a number of prominent sinusoids contained in a 
received signal, estimate the complex amplitude and 
frequency of each of the prominent sinusoids contained in 

10 the received signal. Such methods have applications in the 
field of speech-processing, for example speech recognition 
and speech decomposition. By decomposing a speech signal 
into a set of amplitudes, frequencies and phases of the 
dominant sinusoidal components, speech information can be 

15 efficiently stored or transmitted and then reconstructed 
after being read from a memory, or after being transmitted 
and received. 

Given the already complex nature of this problem 
it is usually assumed that the environment (i.e. channel) 
2 0 and the parameters to be estimated are reasonably 
stationary. In other words, it is assumed that the 
statistical characteristics of the channel and the received 
signal do not change significantly over a short duration of 
time. 

25 Different approaches have been used to address the 

problem of estimating both the complex amplitude and 
frequency of each of a number of prominent frequency 
components contained in a. received signal. The many 
approaches employed include the use of the Discrete Fourier 

30 Transform (DFT) , Adaptive Line Enhancement (ALE) , Extended 
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Kalman Filter (EKF) , Maximum Likelihood Estimation (MLE) , 
and Joint Time Frequency Analysis (JTFA) . 

The DFT was one of the first methods employed 
because it readily enabled the separation of the prominent 
5 sinusoidal components in the frequency domain. This method 
is particularly suitable in applications where the 
parameters are constant, as the DFT can be used in concert 
with the MLE approach. The combination of the DFT and MLE 
approaches results in a method that is the equivalent of 
10 maximizing the periodgram spectrum. Using DFT and MLE 

approaches together, the periodgram can be calculated and 
maximized at discrete frequency points. The prominent 
components of the received signal can be operated upon 
independently, avoiding cross -interference between them. 

15 The combined DFT and MLE method outlined above is 

not very attractive for real time applications because of 
the high computational cost and complexity associated with 
the method. Also, because the method has no memory the 
tracking is not efficient. Several methods have been 

2 0 suggested to enhance the performance of this approach, 

however they have been employed with diminishing returns and 
do not completely resolve the problems. For instance, the 
hidden Markov model was proposed to improve the efficiency 
of tracking the parameters and the computational cost of 
25 this approach can be reduced by use of the Fast Fourier 
Transform (FFT) , in place of the DFT. 

The MLE method on its own is a powerful approach 
to parameter estimation and is widely used in signal 
processing. An iterative method for the MLE extraction of 

3 0 parameters of a harmonic series using the Expectation 

Maximization (EM) algorithm, where the number of harmonics 
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is assumed to be known, has been suggested for AWGN 
(Additive White Gaussian Noise) channels with known noise 
variance. Such methods suffer from poor performance because 
the effect of the cross- interference of harmonics is 
5 ignored. These cross-interferences become significant 

contributors to the degradation of a signal within a low SNR 
(signal-to-noise ratio) environment or for a short data 
length. 

If the characteristic parameters of the signal are 
10 slightly non- stationary or if the additive noise has a time 
varying variance, employing an approach that relies on 
finite length window, such as the DFT, will inherently lead 
to the loss of some information in regard to the dynamics of 
the signal. In previous works using the harmonic model for 
15 a speech signal, the amplitudes, phases and frequencies that 
make-up the speech signal have been estimated. Furthermore, 
it has also been shown that information about the nature of 
the slowly time varying parameters can be obtained from the 
distortions caused by windowing. 

2 0 In yet another approach called Adaptive Line 

Enhancement (ALE) , the frequency of the dominant component 
is estimated by minimizing the output energy of a notch 
filter. An Adaptive Comb Filter (ACF) that is an extension 
of ALE has been used to estimate the frequency of harmonic 

25 signals. 

Kalman Filtering has also been employed in prior 
work for different scenarios. Specifically, EKF has been 
used to track the frequencies, the amplitudes and the phases 
of harmonic components within a periodic signal corrupted by 
30 AWGN. Using similar principles, several non-linear filters 
have been proposed for the decomposition of signals that are 
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modeled as a sum of jointly modulated amplitude and 
frequency cosines in an additive noise environment, where 
the centre frequencies are very slowly time varying. 
Furthermore, assuming that the noise statistics and the 
5 number of superimposed signals are known, an EKF can be 
designed to track frequency formats of speech. 

In another existing approach the signal to be 
analyzed is assumed to be a Polynomial Phase Signal (PPS) 
and unknown parameters are estimated. Several techniques 
10 could be employed to resolve this problem, such as FFT. 
High-resolution frequency estimation methods such as 
Kumaresan-Tuf ts, MUSIC and Matrix Pencil are alternatives to 
estimate the polynomial phase coefficients. 

Exponentially-damped Polynomial Phase Signals 
15 (EPPS) have been treated as a special case of PPS's. In 
such a case, it is typically assumed that a PPS can be 
modeled as having constant amplitude, thus allowing the use 
of JTFA tools such as Wigner-Ville distribution and an 
associated Ambiguity Function. Using the Wigner-Ville 
20 distribution an estimation method that selects an optimal 
time domain window length to resolve the trade-off between 
the estimation bias and the variance of the unknown 
frequency can be used. 

In another existing approach researchers have 
25 considered non- stationary signals as time -dependent ARMA 

(auto-regressive moving average) processes, and suggested a 
general estimation procedure for the ARMA parameters using a 
set of basis functions. In one such work, estimates for the 
signal parameters for any non-stationary signal can be 
3 0 obtained using the Non- linear Instantaneous Least Squares 
(NILS) method. NILS has relation with ML as well as with 
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signal -subspace fitting and linear prediction based 
estimation approaches. 



Summary of the Invention 

According to one broad aspect, the invention 
provides a method of tracking amplitude, phase and frequency 
of a plurality of sinusoidal components in a signal, the 
method comprising: a) processing the signal to produce a new 
set of amplitude and phase estimates using a weighted 
likelihood method; and b) processing the signal to produce a 
new set of frequency estimates using a weighted likelihood 
method. 

In some embodiments, the method further comprises 
sampling the signal to produce a sequence of real -valued 
samples, wherein steps a) and b) are performed in the 
digital domain. 

In some embodiments, the method further comprises 
sampling the signal to produce a sequence of complex-valued 
samples, wherein steps a) and b) are performed in the 
digital domain. 

In some embodiments, steps a) and b) are performed 
in the continuous time domain. 

According to one broad aspect, the invention 
provides a method of tracking amplitude, phase and frequency 
of a plurality of sinusoidal components in a signal, the 
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method comprising: for a current update period: i) 
processing the signal to produce a new set of complex 
amplitude estimates by: a) for a first input set of 
estimated complex sinusoidal components, separating 
5 components to produce component estimates; ii) processing 
the signal to produce a new set of estimated complex 
sinusoidal components by: b) for each component of a second 
input set of estimated complex sinusoidal components, 
estimating a frequency deviation estimate; c) adapting a 
10 previous set of frequency estimates taking into account an 
input set of component estimates and the frequency deviation 
estimates to produce a new set of frequency estimates; and 
d) converting the new set of frequency estimates to a new 
set of estimated complex sinusoidal components. 

15 In some embodiments, the signal is a sequence of 

samples and processing is done in the digital domain. 

In some embodiments, the processing is done in the 
continuous time domain. 

In some embodiments, the further comprises: 
20 performing cross -interference cancellation on the component 
estimates to produce a new set of cross- interference 
cancelled component estimates, and using the new set of 
cross-interference cancelled estimates as the input set of 
component estimates in an execution of step c) . 

25 In some embodiments, the method further comprises: 

performing complex envelope extraction on the component 
estimates to produce a new set of complex amplitude 
estimates . 

In some embodiments, the method further comprises: 
30 performing complex envelope extraction on the cross- 
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interference cancelled component estimates to produce a new 
set of complex amplitude estimates. 

In some embodiments, for the first input set of 
estimated complex sinusoidal components, separating 
components to produce component estimates is done using a 
weighted log-likelihood function with a first weighting 
sequence; for each of the second input set of estimated 
complex sinusoidal components, estimating the frequency 
deviation estimate is done using a weighted log-likelihood 
function with a second weighting sequence. 

In some embodiments, the first and second 
weighting sequences are the same. 

In some embodiments, step i) comprises a first 
half -iteration, and step ii) comprises a second half 
iteration, one first half -iteration and one second half- 
iteration comprising a complete iteration and wherein for 
each update period, a plurality of complete iterations are 
performed to produce the new set of complex amplitude 
estimates and the new set of estimated complex sinusoidal 
components . 

In some embodiments, the first input set of 
estimated complex sinusoidal components and the second set 
of estimated complex sinusoidal components are initially set 
to initial values, and thereafter are set to estimated 
25 complex sinusoidal components produced by a previous 
iteration of the method. 

In some embodiments, for each update of the 
complex amplitude and frequency: the step of processing 
samples of the sequence of samples to produce a new set of 
30 complex amplitude estimates is performed before the step of 
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processing the sequence of samples to produce a new set of 
estimated complex sinusoidal components; the first input set 
and the second input set of estimated complex sinusoidal 
components comprise the new set of estimated complex 
5 sinusoidal components determined during a previous update 
period; wherein the input set of cross- interference 
cancelled estimates comprises the new set of cross- 
interference cancelled estimates determined during the 
current update period. 

10 In some embodiments, for each update of the 

amplitude, phase and frequency: the step of processing the 
signal to produce a new set of estimated complex sinusoidal 
components is performed before the step of processing the 
sequence of samples to produce a new set of complex 

15 amplitude estimates; the input set of component estimates 
comprises the set of cross-interference cancelled estimates 
determined during a previous update period; the first input 
set of estimated complex sinusoidal components comprises the 
new set of estimated complex sinusoidal components 

2 0 determined during the current update period and the second 
input set of estimated complex sinusoidal components 
comprises the new set of estimated complex sinusoidal 
components determined during a previous update period. 

In some embodiments, for the first input set of 
25 estimated complex sinusoidal components, performing 
component extraction using a weighted log-likelihood 
function with the first weighting sequence comprises 
filtering the samples with a respective component extraction 
filter tuned to a respective one of the first input set of 
30 estimated complex sinusoidal components. 
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In some embodiments, performing cross-interference 
cancellation on the component estimates to produce a new set 
of cross- interference cancelled component estimates 
comprises multiplying the component estimates by a cross- 
5 interference cancellation matrix. 

In some embodiments, performing complex envelope 
extraction on the cross-interference cancelled component 
estimates to produce the new set of complex amplitude 
estimates comprises multiplying each cross-interference 
10 cancelled component estimate by the respective estimated 
complex sinusoidal component with negative exponent. 

In some embodiments, for each of the second input 
set of estimated complex sinusoidal components, estimating a 
frequency deviation estimate using the weighted log- 
15 likelihood function with the second weighting sequence 

comprises filtering the sampled sequence with a respective 
frequency deviation filter tuned to the estimated complex 
sinusoidal component . 

In some embodiments, adapting the previous set of 
20 frequency estimates taking into account an input set of 

component estimates and the frequency deviation estimates to 
produce a new set of frequency estimates comprises applying 
an adaptation value to each previous frequency estimate, the 
adaptation value being a function of both the input set of 
25 component estimates and the frequency deviation estimates. 

In some embodiments, applying an adaptation value 
to each previous frequency estimate, the adaptation value 
being a function of both the input set of component 
estimates and the frequency deviation estimates comprises: 
30 determining a partial derivative with respect to each 

^ r 
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estimated complex sinusoidal component of a function based 
on the weighted log-likelihood function; for each frequency 
estimate, determining the adaptation value from the 
respective partial derivative. 

5 In some embodiment, adapting the previous set of 

frequency estimates taking into account the input set of 
component estimates and the frequency deviation estimates to 
produce the new set of frequency estimates comprises: 
applying an adaptation value to each frequency estimate in 

10 the previous set of frequency estimates, the adaptation 

value being a function of both the component estimates and 
the frequency deviation estimates to produce an intermediate 
set of frequency estimates; using the frequency deviation 
estimates and previous frequency deviation estimates to 

15 produce an estimate of chirp for each sinusoidal component; 
for each sinusoidal component, combining the frequency 
deviation estimate and the estimate of chirp to produce a 
new frequency estimate. 

In some embodiments, converting the new set of 
20 frequency estimates to new estimated complex sinusoidal 
components comprises combining previous estimated complex 
sinusoidal component estimates with the new frequency 
estimates . 

In some embodiments, combining the previous 
25 estimated complex sinusoidal component estimates with the 

new frequency estimates comprises: multiplying each previous 
estimated complex sinusoidal component estimate by e A (j x 
new frequency estimate) . 
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In some embodiments, one or more ASICs 
(application specific integrated circuit) adapts to 
implement a method. 

In some embodiments, one or more DSPs (digital 
signal processors) adapts to implement a method. 

In some embodiments, one or more FPGAs (field 
programmable gate arrays) adapts to implement a method. 

In some embodiments, one or more general purpose 
processors adapts to implement a method. 

In some embodiments, a combination of at least two 
circuits selected from a group consisting of ASIC, FPGA, 
DSP, and general purpose processor adapts to implement a 
method. 

In some embodiments, a computer readable medium 
having executable code embodied therein for causing a 
processing platform to execute a method. 

According to one broad aspect, the invention 
provides an apparatus for tracking amplitude, phase and 
frequency of a plurality of sinusoidal components in a 
signal, the apparatus comprising: a first processing path 
adapted to process the signal to produce a new set of 
amplitude and phase estimates using a weighted likelihood 
method; and a second processing path adapted to process the 
signal to produce a new set of frequency estimates using a 
weighted likelihood method. 

In some embodiments, the apparatus further 
comprises: a sampler adapted to sample the signal to produce 
a sequence of real -valued samples, wherein the first and 
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second processing paths perform signal processing in the 
digital domain. 

In some embodiments, an apparatus further 
comprises: a sampler adapted to sample the signal to produce 
5 a sequence of complex- valued samples, wherein the first and 
second processing paths perform signal processing in the 
digital domain. 

In some embodiments, the first and second 
processing paths perform signal processing in the continuous 
10 time domain. 

According to one broad aspect, the invention 
provides an apparatus for tracking amplitude, phase and 
frequency of a plurality of sinusoidal components in a 
signal, the apparatus comprising: at least one component 

15 extraction filter adapted process the signal to produce 
component estimates for each of a first input set of 
estimated complex sinusoidal components, each component 
extraction filter being tuned to a respective one of the 
first input set of estimated complex sinusoidal components; 

20 at least one frequency deviation filter adapted to process 
the signal to produce a frequency deviation estimate for 
each of a second input set of estimated complex sinusoidal 
components, each frequency deviation filter being tuned to a 
respective one of the second input set of estimated complex 

25 sinusoidal components; at least one adaptive frequency 

tracker adapted to produce a new set of frequency estimates 
by adapting a previous set of frequency estimates taking 
into account an input set of component estimates and the 
frequency deviation estimates; and at least one component 

30 generator adapted convert the new set of frequency estimates 
to a new set of estimated complex sinusoidal components. 
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In some embodiments, the signal is a sequence of 
samples and processing is done in the digital domain, and 
wherein the at least one component generator comprises at 
least one digital controlled oscillator. 

5 In some embodiments, the apparatus further 

comprises: a cross-interference canceller adapted to perform 
cross-interference cancellation on the component estimates 
to produce a new set of cross-interference cancelled 
component estimates; wherein the new set of cross- 
10 interference cancelled estimates are used as the input set 
of component estimates to the adaptive frequency tracker. 

In some embodiments, the apparatus further 
comprises: at least one complex envelope estimator adapted 
to perform complex envelope extraction on the component 
15 estimates to produce a new set of complex amplitude 
estimates . 

In some embodiments, the apparatus further 
comprises: at least one complex envelope estimator adapted 
to perform complex envelope extraction on the cross- 
20 interference cancelled component estimates to produce a new 
set of complex amplitude estimates. 

In some embodiments, each component extraction 
filter implements a weighted log-likelihood function with a 
first weighting sequence; each frequency deviation filter 
25 implements a weighted log- likelihood function with a second 
weighting sequence. 

In some embodiments, the first and second 
weighting sequences are the same. 
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In some embodiments, the first input set of 
estimated complex sinusoidal components and the second set 
of estimated complex sinusoidal components are initially set 
to initial values, and thereafter are set to previously 
5 determined estimated complex sinusoidal components. 

In some embodiments, for each time a new set of 
complex amplitude estimates is produced by the apparatus: 
the component extraction filter (s) operate to produce the 
new set of complex amplitude estimates before the frequency 

10 deviation filter (s) operate to produce the new set of 

estimated complex sinusoidal components; the first input set 
and the second input set of estimated complex sinusoidal 
components comprise the new set of estimated complex 
sinusoidal components determined during a previous update 

15 period; wherein the input set of cross- interference 
cancelled estimates comprises the new set of cross- 
interference cancelled estimates determined during the 
current update period. 

In some embodiments, for each time a new set of 
20 complex amplitude estimates is produced by the apparatus: 
the component extraction filter (s) operate to produce the 
new set of estimated complex sinusoidal components before 
the frequency deviation filters operate to produce the new 
set of complex amplitude estimates; the input set of 
25 component estimates comprises the set of cross-interference 
cancelled estimates determined during a previous update 
period; the first input set of estimated complex sinusoidal 
components comprises the new set of estimated complex 
sinusoidal components determined during the current update 
3 0 period and the second input set of estimated complex 

sinusoidal components comprises the new set of estimated 
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complex sinusoidal components determined during a previous 
update period. 

In some embodiments, the cross- interference 
canceller produces the new set of cross-interference 
5 cancelled component estimates by multiplying the component 
estimates by a cross-interference cancellation matrix. 

In some embodiments, the complex envelope 
estimator (s) produce the new set of complex amplitude 
estimates by multiplying each cross-interference cancelled 
10 component estimate by the respective estimated complex 
sinusoidal component with negative exponent. 

In some embodiments, the adaptive frequency 
tracker (s) apply an adaptation value to each previous 
frequency estimate, the adaptation value being a function of 
15 both the component estimates and the frequency deviation 
estimates . 

In some embodiments, the adaptive frequency 
tracker (s) determine a partial derivative with respect to 
each estimated complex sinusoidal component of a function 
20 based on a weighted log-likelihood function and for each 

frequency estimate, determine the adaptation value from the 
respective partial derivative. 

In some embodiments, the adaptive frequency 
tracker (s) produce a new set of frequency estimates by 

25 applying an adaptation value to each frequency estimate in a 
previous set of frequency estimates, the adaptation value 
being a function of both the component estimates and the 
frequency deviation estimates to produce an intermediate set 
of frequency estimates, and using the frequency deviation 

30 estimates and previous frequency deviation estimates to 
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produce an estimate of chirp for each sinusoidal component, 
and for each sinusoidal component combine' the frequency 
deviation estimate and the estimate of chirp to produce a 
new frequency estimate. 

In some embodiments, the component generator (s) 
convert the new set of frequency estimates to new estimated 
complex sinusoidal components by combining previous 
estimated complex sinusoidal component estimates with the 
new frequency estimates. 

In some embodiments, a computer in combination 
with a computer readable medium compatible with the computer 
are provided which are cooperatively adapted to implement 
any of the above methods . 

Brief Description of the Drawings 

Preferred embodiments of the invention will now be 
described in greater detail with reference to the 
accompanying diagrams, in which: 

Figure 1 is a block diagram of an apparatus for 
the adaptive estimation and tracking of a multi -component 
sinusoidal real -valued observed signal provided by an 
embodiment of the invention; 

Figure 2 is a flow chart of a method provided by 
an embodiment of the invention for the adaptive estimation 
and tracking of a multi -component sinusoidal real -valued 
observed signal; 

Figure 3A is a block diagram of a first Component 
Extraction Filter (CEF) usable in the apparatus of Figure 1; 



73674-4 



18 

Figure 3B is a block diagram of a second CEF 
usable in the apparatus of Figure 1; 

Figure 4A is a block diagram of a first Frequency 
Deviation Filter (FDF) usable in the apparatus of Figure 1; 

5 Figure 4B is a block diagram of a second FDF 

usable in the apparatus of Figure 1; 

Figure 5 is a block diagram showing details of the 
CIC and CEEs of Figure 1; 

Figure 6 is a flow chart of a method for 
10 dynamically detecting and updating the number of prominent 
sinusoidal components to be tracked in the real -valued 
observed signal; 

Figure 7 is a block diagram of a speech coder 
employing the adaptive estimation and tracking method of 
15 Figure 2; 

Figure 8 is an example implementation of one of 
the DCOs of Figure 1; 

Figures 9A and 9B contain plots of frequency 
tracking for different SNRs for a signal with two 
20 components. Estimated frequencies (solid) and true 

frequencies (dotted) are superimposed on background of the 
spectrum of the main signal, 9A: Tracked frequencies of the 
first scenario for SNR=0dB, 9B: SNR=20dB; 

Figure 10 contains plots of estimation errors of 
25 the proposed algorithm; solid: mean squared error \x n -x n \ 2 
(averaged over 50 runs), dotted: sum of MSE of components 
I - A I 2 I _ * I 2 
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Figure 11 contains plots of estimated amplitudes 
in 20dB, using a Hamming window with a length of 129, 
dotted: true values, solid and dashed: estimated values; 

Figure 12 contains plots of average of squared 
5 frequency estimation error of the proposed algorithm: 

2 

• (averaged over 50 runs) . Solid: one 



A-A 2 + f 2 -f 2 



2 

iteration per time instance with p. = 1CT 4 , Dashed: two 
iterations per time instance ^ = 10" 4 and \i 2 = 5 x 10" 5 ; and 

Figures 13A and 13B contain plots of results of 
10 decomposing a segment of speech to four components using a 
47 sample length Hamming window, two iteration and different 
y for the components, 13A: Tracked frequencies on background 
of the spectrogram of the main speech, 13B: Spectrogram of 
the constructed signal. 

15 Detailed Description of the Invention 

The problem to be solved can be conceptualized in 
a discrete time mathematical model. In such a model, the 
signal to be examined can be considered a real -valued 
observed signal x n having L sinusoidal AM- FM components and 

20 corrupted by additive white noise. The real-valued observed 
signal x n consists of a sequence of real -valued samples of a 
multi -component real-valued signal. The samples are 
obtained by sampling with a sampling period of Tand sampling 
frequency f s = 1/7 Hz. The real-valued observed signal can be 

25 represented mathematically as follows: 

x.^R&ia^ + N.; neZ (1) 
/=i 
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where N n e 91 is a real additive white noise sample, the 
complex signal a t = \a\e J(Pl eC represents the amplitude and phase 
of the / th component to be estimated, and 0) { G[-7t,7r] is the 
frequency (radian/sample) of the / th component to be 
5 estimated. Re ( . ) denotes the real part of ' a complex number. 
In equation (1) , it is noted that when changing the pair 
{co n aj) to (-a) n a { ), where (.)* is the complex conjugate 
operator, the observed signal is unchanged. Thus, the sign 
of the frequency is not identifiable. Because of this, it 
10 can be assumed that <o l e[0,^] , and if the method used for 
estimation results in a negative value for the radian 
frequency o l the pair {co n a } ) can be changed to (-0) n a*) without 
a loss of any information. The relationship between 0) l and 

the real continuous frequency f s is given by /J=— /. 

2n 

15 The problem that is addressed here is to estimate 

the complex amplitude (i.e. amplitude and phase) and the 
frequency of all prominent sinusoidal components of such a 
signal . 

It is noted that an embodiment of the invention 
20 also provides a similar but simpler solution for the case of 

complex observed signals, where ^=^ =I fl/ ¥ +^6C. The 

solution for the complex case is a simplified version of the 
real case solution, since in the complex case the quadrature 
components of the signal are also observed. 

25 It is assumed that the amplitudes, the phases and 

the frequencies of the prominent sinusoidal components are 
very slowly time-varying or equivalently they are assumed to 
be band limited and smooth signals. Furthermore, it is 
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assumed that each of the prominent sinusoidal components may 
disappear or appear, but does so in such a manner that the 
number of prominent sinusoidal components rarely changes. 
It is also assumed that the number of prominent sinusoidal 
5 components is known initially. The method may still work 
should one or more of these assumptions fail. 

A likelihood function can be evaluated which 
represents the amount of information about the received 
(observed) signal that is available to the receiver. 
10 Evaluation of the likelihood function can provide an 

estimation of the unknown parameters. If it is assumed that 
N n is a zero-mean white Gaussian random process with 
variance cr 2 N (n) , the log-likelihood function at time n of the 
observed x n can be expressed as follows: 



the maximization of the above likelihood function over a 
rectangular window (e.g. by expectation maximization 

20 method) could provide an appropriate estimate for the 

characteristic parameters of each of the prominent sinusoids 
contained in the observed real-valued signal. However, this 
approach requires a large amount of computation and also 
suffers from the problem of local optima which can result in 

25 inaccurate results. 



15 U 




It is noted that within a stationary environment 



Embodiments of the invention provide a method and 
system for evaluating the likelihood function that is 
computationally feasible and provides accurate estimates for 
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the characteristic parameters of each of the prominent 
sinusoids contained in the observed real-valued signal. 



Referring now to Figure 1, shown is a block 



diagram of an apparatus provided by an embodiment of the 
5 invention for the adaptive estimation and tracking of a 

multi -component real-valued observed signal x n . It is noted 
that in what follows, two different indices for time are 
used. The index u n" is used to refer to the processing 
performed by the method at current time n. The index w i" is 
10 a dummy variable used to refer to the observed signal at 
times other than the current time n. Typically, the 
processing at time n uses multiple different observed 
signals 



15 is a set of Component Extraction Filters (CEFs) 110. There 
is one CEF for each prominent sinusoidal component being 
tracked. It is assumed there are L such components where 
L>1. The CEFs 110 have a first input 111 and a second input 
112. The first input 111 accepts the real-valued observed 

2 0 signal x n , while the second input 112 accepts a set of 
estimated complex sinusoidal components 



P")/ = 1,. J Z corresponding to the prominent frequency 

components in the real -valued observed signal x\ being 
tracked. The CEFs 110 also have an output 113 from which 

25 they deliver a set of component estimates Y n made up of 
estimates of the prominent signal components of the real- 
valued observed signal x n . Contained within CEFs 110 are a 
number of filters, preferably band pass filters, each of 
which is tuned to a respective frequency corresponding to 

3 0 one of the prominent signal components of the real -valued 



A first block of the apparatus, shown in Figure 1, 
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observed signal x n . More specifically, in one embodiment at 
time n these are L band pass filters having impulse 
responses {/*,„}/ = 1, 2, ..L . These band pass filters are 
described in further detail below. 

5 A second block of the apparatus is a Cross- 

Interference Canceller (CIC) 130. The CIC 130 has a first 
input 131 and a second input 132. The first input 131 
accepts the set of component estimates Y n , while the second 
input 132 accepts the set of estimated complex sinusoidal 
10 components The CIC 13 0 also has an output 13 3 from 

which it delivers a set of cross- interference cancelled 
component estimates X n of the prominent signal components of 
the real -valued observed signal x n . 

A third block of the apparatus is a set of Complex 
15 Envelope Estimators (CEEs) 150. The CEEs 150 have a first 
input 151 and a second input 152. The first input 151 
accepts the set of cross- interference cancelled signal 
component estimates X n , while the second input 152 accepts 
the set of estimated complex sinusoidal components The 
20 CEEs 150 also have an output 153 from which they deliver a 
set of complex amplitude estimates {a ln }. The CEEs 150 are 
made up of a number of signal multipliers, each of which is 
usable for a respective frequency corresponding to one of 
the prominent signal components of the real -valued observed 
25 signal x n . 

A fourth block of the apparatus is a set of 
Frequency Deviation Filters (FDFs) 120. The FDFs 120 have a 
first input 121 and a second input 122. The first input 121 
accepts the real-valued observed signal x n , while the second 
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input 122 accepts the set of estimated complex sinusoidal 
components {p"}. The FDFs 120 also have an output 123 from 

which they deliver a set of frequency deviation estimates Y n 
each a measure of a frequency deviation of a prominent 
5 signal component of the real -valued observed signal x n . The 
FDFs 120 consist of a set of L filters, each of which is 
tuned to a unique frequency corresponding to one of the 
prominent signal components of the real -valued observed 
signal x n . The filters have impulse responses {/*, „}/ = 1, 2, ...,Z . 

10 A fifth block of the apparatus is a set of Digital 

Controlled Oscillators (DCOs) 140. The DCOs 140 have an 
input 141 and an output 142. The input 141 accepts a set of 

frequency estimates ^ corresponding to the frequencies 
of prominent frequency components contained in up the real- 
15 valued observed signal x n . The output of the DCOs 140 
available at the output 142 is a new set of estimated 

complex sinusoidal components jf^} produced by combining the 
previous estimated complex sinusoidal components with 
the frequency estimates 

20 A sixth and last block of the apparatus that is 

shown in Figure 1 is a set of Adaptive Frequency Trackers 
(AFTs) 160. The AFTs 160 have a first input 161 and a 
second input 162. The first input 161 accepts the set of 

frequency deviation estimates Y n , while the second input 162 
25 accepts the set of cross-interference cancelled component 

estimates X n . The AFTs 160 also have an output 163 from 

which they deliver the set of frequency estimates {<££,„}. 
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It is noted for the single component case, only 
one CEF, CEE, FDF, DCO ad AFT are required and there is a 
single CIC so long as there are two or more components, no 
CIC being required for the single component case. 

5 In operation, the real-valued observed signal x n is 

simultaneously fed to the CEFs 110 and the FDFs 120 via 
inputs 111 and 121 respectively. The adaptive joint 
estimation and tracking method provided by the invention is 
recursive, and the method may equivalently start at either 

10 the CEFs 110 or the FDFs 120. In alternate half -iterations 
of the method, either the complex amplitude estimates {a ln } 
are updated as a function of the observed signal x n and 
previous estimated complex sinusoidal components {q"|, or the 
estimated complex sinusoidal components {q"| (and the 

15 frequency estimates [cb Kn )) are updated as a function of the 
observed signal x n and previous cross- interference cancelled 
estimates X n . The CEFs 110 and the FDFs 120 generate values 
for the component estimates Y n and frequency deviation 

estimates Y n respectively. Both the CEFs 110 and the FDFs 
20 120 are initialized with estimates of the estimated complex 
sinusoidal components, corresponding to each of the 
prominent sinusoidal components contained in the real -valued 
observed signal x n . Furthermore, both the CEFs 110 and the 
FDFs 120 indirectly supply the DCOs 140 a feedback signal 
25 that allows the DCOs 140 to update the frequency estimates. 

The CEFs 110, CIC 130 and CEEs 150 collectively 
process the observed signal x n to produce complex amplitude 
estimates a ln of the prominent sinusoidal components. Each 
of the filters in CEFs 110 filters the observed signal x n to 
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produce a respective initial estimate Y n of each frequency 
component. The filters are tuned to look at frequencies 
specified by the estimated complex sinusoidal components 

{d"}. The CIC 130 accepts the component estimates Y n 
5 generated by the CEFs 110 and the estimated complex 

sinusoidal components generated by the DCOs 140. The CIC 
130 can be basically described as a matrix processor that 
combines the estimated complex sinusoidal components with 
the component estimates Y n to produce the cross -interference 

10 cancelled component estimates X n . The mathematical details 
of this block will be given in detail in what follows. 

The CEEs 150 operate by multiplying the estimated 
complex sinusoidal components and the corresponding cross- 
interference cancelled component estimates X n to produce the 
15 set of complex amplitude estimates, each complex amplitude 
estimate corresponding to a respective prominent sinusoidal 
component contained in the observed signal x n . 

The FDFs 120, DCOs 140 and AFTs 160 collectively 
process the observed signal x n to produce the estimated 

20 complex sinusoidal component jfi"} and frequency estimates 
The FDFs 120 generate estimates of the prominent 
sinusoidal components frequency deviations from the real- 
valued observed signal and previous estimates of the 
frequencies present in the real -valued observed signal x n . 

25 The filters in the FDFs 120 filter the observed signal to 
produce the set of estimates Y n of deviations in the 
modulated frequency of the prominent sinusoidal components 
from the previous estimates. The FDFs 120 pass the 

frequency deviation estimates Y n they have generated to the 
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AFTs 160 as shown in Figure 1. The AFTs 160 use the 
frequency deviation estimates Y n from the FDFs 120 and the 
cross -interference cancelled component estimates X n from the 
CIC 130 to generate the set of frequencies {o) ln } 

5 corresponding to the respective components contained the 
real -valued observed signal. 

The frequencies \S ln } generated by the AFTs 160 are 
output by the method as the new set of frequency estimates. 
They are also passed to the DCOs 140, which contain 

10 recursive complex sinusoidal signal generators that modulate 
the frequency estimates by combining them with previous 
estimates to produce estimated complex sinusoidal components 
and in so doing reduces the effect of noise in the frequency 
estimates. Filtering the frequency estimates in this way 

15 minimizes the amount of computational error propagated by 

erroneous frequency estimation. As previously indicated the 
DCOs 14 0 then send the new estimated complex sinusoidal 
components {ty} to the CEFs 110, the FDFs 120, the CIC 130 

and the CEEs 150, so that the next iteration of processing 
2 0 can proceed. 

The mathematics behind the approach embodied in 
Figure 1 will now be presented. The adaptive estimation and 
tracking method provided by an embodiment of the invention 
is obtained by maximizing a weighted average of equation 
25 (2), the likelihood equation given above. The weighted- 
average of equation (2) is taken over a finite amount of 
time, or rather within a window in time. A recursive 
adaptive filter with a low order of computational complexity 
is employed to achieve this. 
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By appropriately selecting the averaging window 
the problem of converging on an erroneous local optima (e.g. 
what is known as a false- lock in a PLL) can be overcome, 
with the additional benefit of controlling the noise 
5 limiting and tracking behavior of the adaptive estimation 
and tracking method. In order to estimate the parameters 
characterizing prominent sinusoidal components comprising 
the real-valued observed signal x n at a time n, in a 
preferred embodiment the following parametric windowed 
10 (weighted) log-likelihood function is maximized: 

L{n) = £ w ^ifc(0.®,(0.^(0L) < 3 > 

J=-CO 

In equation (3) W k is a window function described in detail 
below, where k is an index for the window function which is 
set relative to n, i. 

15 This approach will be referred to herein as the 

Maximum Weighted Likelihood (MWL) approach. Despite being 
presented within a system that assumes a stationary 
environment, it should be noted that this approach is also 
suitable for non- stationary and adaptive parameter 

20 estimation. 

Furthermore, the weighted log-likelihood function 
can be employed, and adapted for implementation using the 
apparatus of Figure 1, by defining the impulse response \h ln ) 

and ]/*/,„} as described below. In another embodiment, the 
25 structure of Figure 1 is provided with no specific 
constraints on the impulse responses {h ln } and {ft/,*}. 

Referring to equation (3), the value W n _. is the 
weight of the information received at time n-i in order to 
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estimate the unknown parameters at time n. In some 
embodiments, the window function, \N k , is selected to satisfy 
one or more of the following conditions: 

1. The window is normalized (i.e. Z^W^=1). 

5 2. The window is centered at zero (i.e. SkW^=0) . 

3. The causality of the adaptive estimation and 
tracking method relies on > D => = 0, where D 
denotes an acceptable delay for the estimated 
parameters characterizing the prominent sinusoidal 

10 components comprising the real -valued observed 

signal x n . 

4. The condition for L(n) to yield a positive 
integration from the information gathered at 
different time instances is that: W k > 0 . 

15 The first condition is not absolutely necessary. 

It has been included because it simplifies the discussion in 
regard to the adaptive estimation and tracking method being 
described. The second condition could also be violated, in 

such a case, the value of X^L,^* represents the delay lag 

20 of the estimation. Again the second condition has been 

included because it simplifies the upcoming description of 
the adaptive estimation and tracking method. More 
generally, if one or more of the conditions are not 
satisfied, the performance might degrade depending on the 

25 situation. 



Considering more carefully the above conditions, 
it can be observed that if is considered to be the unit 
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impulse response of a linear time -invariant filter, this 
filter is strictly a lowpass filter having its maximum 
frequency domain gain for the DC component. 



Additionally, if the window weights are all 



5 non-negative, L(n) represents a measure of the received 

information about the parameters in the neighborhood of time 
n. This window will determine the resolution of the 
frequency estimation and the extracting filters for the 
complex amplitude a { . 

10 The maximization of the criterion in equation (3) 

does not resolve the problem because the number of unknown 
parameters in equation (3) is more than the number of 
parameters obtained by observation within the window. So as 
an approximation, it is assumed that the unknown parameters 

15 are approximately constant across the course of the window 
and this approximation allows the likelihood function to be 
simplified. Consequently, the following simplified 
likelihood function is maximized in place of equation (3) : 



information about the real-valued observed signal. The 
approximation error reduces as the variation of the observed 
25 parameters reduces as caused by shortening the window 

length. Yet with a longer window length a longer duration 
of the observed signal can be used. This has the effect of 





(4) 



20 



Due to the approximation made above an inherent 
approximation error is introduced resulting from ignoring 



73674-4 



31 

increasing the estimation accuracy as long as the variance 
of the observed parameters is sufficiently small. However, 
it is not always the case that the variances of the 
parameters characterizing the prominent sinusoidal 
5 components of the real-valued observed signal will be 
sufficiently small and such circumstances can not be 
guaranteed to occur in practical applications of the 
invention. This shows a trade-off to be considered when 
choosing the duration of the time domain window. The length 

10 of the time domain window is preferably chosen so that the 
effects of the varying parameters characterizing the 
observed (received) signal are balanced against the 
requirement to observe an adequate portion of the observed 
signal. Thus, it is clear that a shorter time domain window 

15 length is preferable in practical embodiments of the 

invention. The exact length of the time domain window may 
be determined for a particular application by empirical 
methods . 

Maximizing L{n) with respect to cr 2 N (n) , the result 
20 is the following equation that approximates the variance: 



= -Et, Re («/e M f W„, (5) 



Next, the value of d 2 N (n) is substituted from 

equation (5) into equation (4) and then L is maximized with 
25 respect to other parameters. The result of this operation 
is another estimate of the likelihood function as follows: 



(6) 
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Maximizing equation (6) yields to minimization of 
the following Weighted-Least-Square (WLS) criterion for 
estimation of amplitude and frequency at time m 



/=-oo 



2 

W. (7) 



5 The window has a direct impact on the lock- in 

range of the adaptive estimation and tracking method, 
wherein the main-lobe width of the Fourier spectrum of 
defines the lock- in range. In general the lock- in range of 
a method is the maximum initial frequency deviation that can 
10 be tolerated by the method such that the method can acquire 
and begin to track the input frequency. 

Throughout the rest of this disclosure, the 
following notations and definitions are adopted in order to 
provide all of the mathematical details of the present 
15 invention: 

the vector of the complex modulated amplitudes: 



X n ^[a x e^\-,a L e^"] 

= [ X \,n > * ' ' > X L,n ] = n + jX i n 



the estimation of the clean signal: 
20 the estimation error: 



^ n X n X n 



the derivative of a real-valued function f(v) with respect to 
a complex variable v: 
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df{y) _ 1 df(v) j df(y) 



dv 2 dRe(v) 2 dlm(v) 
A window: w ( . > 0 and W(l) = 1 : 



To begin, a detailed derivation of how the complex envelope 
estimates {a ln } are computed by blocks 110, 13 0, 150 of Figure 
1 will be presented. 



10 To maximize J(n) set ^ = 0 for 1 = 1,..., L. Thus, 



= ZZRefa^^^^-'^w,. - Z^e^y (8) 



j=-co p-\ 

= It— ^ -«*" Xv w w ( . 



/=-oo p~\ 

Given the values of the frequencies {tfJ,}^, the 
above system (equation (8)) is solved to provide estimates 
of the complex amplitudes {tf p }^ =1 for each of the prominent 

15 sinusoidal components. To simplify the solution of the 
above system, it is useful to define the following set of 
linear filters: 



73674-4 



34 

IX*"' =W(e^' Z ) = W(Q,z). (9) 

Z=-00 

Component estimate (collectively the vector Y n 

of Figure 1) can be computed according to: 

/=-oo 

5 Referring to Figure 3A and the filters given by 

equation (9) , it can be seen that y ln can be computed as the 
output of a linear bandpass filter 200 that is tuned to the 
estimated centre frequency & l , for the I th prominent 
sinusoidal component of the real-valued observed signal. 

10 The CEFs 110 of Figure 1 include one such band pass filter 
200 for each component. Figure 3B shows an equivalent 
lowpass filter implementation of the filter shown in Figure 
3A. A multiplier 208 multiplies the input x n by e j6),n . The 
result is lowpass filtered with a lowpass filter 210 having 

15 an impulse response equal to the weighting function w k . The 

output is multiplied by e~ j(D,n with multiplier 212. Both 
filters, shown in Figures 3A and 3B, attenuate all 
components except the corresponding component and the 
corresponding component appears at the output with unity 
20 gain. In other words the filter is substantially 

dimensionless for a particular component (// / (Q / ) = l) when the 
estimated frequency is accurate. In this manner, it is easy 
to understand how the adaptive estimation and tracking 
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method operates to reject additive noise. The filter H { (z) 
should be designed to pass through the / th component. The 
bandwidth of H,(z) equals the bandwidth of W and should be 
more than the bandwidth of the amplitudes. 

Each filter is only adjusted by one parameter, 
namely the previous estimated complex sinusoidal component, 
and the output of each filter will contain some interference 
from other components arising from non-zero gain in the stop 
band of each filter. These cross-interferences can be 
removed by summarizing (8) in the following equation: 



P>J(n\ l *> a e m + m ' Xm - i) + a *e Jie¥ ''* H ~ i) 



p=\ i=- 



= e 1<0 ' n 



(11) 



From — — - = 0, the following is derived 
da l 

2*, = Zfl^W^^J+ay^W^"^ (12) 
P =\ 

To have a matrix form solution, it is useful to 
define the following: 

Y n =fou„>-,2y L J =2Y r . n+ 2jY. n 

MoL-wfon,) d3) 

[^)L-w(fi,/fiJ 

Note that there is no need for a computational 
division operation in calculation of W(q / /Q /7 ) as 

w(q,/qJ= \A/(q,Q*J. Thus from equation (12) and the 
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definitions given in (13) , the linear system can be solved 
by solving the following complex set of equations: 

Y n =w(n)x n +v(n)x: a*) 

where (.)* denotes the conjugate operator. Expanding the L- 
5 dimensional complex-valued system of (14) into a 2L- 

dimensional real-valued system, the following is the result: 



X" 









Re{w(n) + K(b)} - \m{w(Q) - v(q)} 
\m{w(n) + v(q)} Re{^(Q) - v{n)} 



(15) 



When the frequencies of the prominent sinusoidal 
components are far enough apart F(q) will not be ill 
10 conditioned. 

Thus, by applying the inverse of F(q) the 
following solution determines the cross- interference 
cancelled component estimates 



X, 



= F Q 



(16) 



15 F _, (q) removes the cross- interference of adjacent 

prominent sinusoidal components and will be referred to as 
the "cross- interference cancelled matrix", or CIC matrix. 
Using (16) , the cross-interference cancelled modulated 
component estimates (i.e. the vector X n =X rn +jjt in ) are 

20 obtained. As J(n) is a quadratic function of the amplitudes 
the above solution is the optimum. 
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The complex amplitudes of the prominent sinusoidal 
components are then estimated by: 

i =[^,a u ,-,6j = [e-^,e-^\-,e-^"] *X„ (17) 

where the symbol * denotes element-wise multiplication of 
the elements of two vectors. Referring now to Figure 1, the 
CIC 13 0 operates to multiply the components output by the 
CEFs 110 by the CIC matrix to produce the cross-interference 

cancelled components X n . The CEEs 150 calculate equation 

(17) and output the complex envelope estimates d t . The 

operation of the CIC 130 and CEEs 150 together is shown in 
further detail in Figure 5, which shows the component 
estimates before CIC y n ={y />#1 } input to a cross-interference 

cancellation matrix multiplication which multiplies Y n by 

F' ] (Cl). The output X n ={x ln } is multiplied by the set of e' s&,n 

in the multipliers 152 of the CEEs 150 to produce the 
complex amplitude estimates. 

The optimization of J(n) to track and to estimate 
the frequencies of the prominent sinusoidal components 
involves non- linear time varying filtering. According to a 
20 preferred embodiment of the invention, an adaptive method 
for the estimation and tracking of frequencies and 
amplitudes is provided. 

It is assumed that an accurate initial value is 
provided for the frequency of each prominent sinusoidal 
25 component. This initial value can be obtained from the 

previous iteration of the method, by another method or by a 

simple initialization procedure. Then, once X n is obtained 



10 



73674-4 



38 

from (16) it is substituted into J(n) as shown in the 
following : 



/=! /=! 



An) = X 



1=-C0 / = 1 



(18) 



Following the above substitution yielding equation 
5 (18) J(n) is minimized with respect to the unknown 

frequencies, each of the unknown frequencies corresponding 
to each of the prominent sinusoidal components. Once this 

is completed, an iterative approach may be used to find X n 
again and so on until the desired level of accuracy is 

10 obtained. However, for this example, only a single 

iteration for each time step n will be demonstrated. The 
adaptive methods work with only one iteration. Increasing 
the number of iterations provides some improvement at the 
expense of increased computational complexity. The 

15 frequency estimates obtained by the adaptive method for the 
estimation and tracking method is then used as the initial 
value (estimate) for the next time iteration n+I. 



To optimize, one may ignore the dependency between 



X n and Q and compute 



am 

da 



as follows 



20 



8J(n) 
da>, 



= -lU-> - ))mfe, fl e--«w j ) 



(19) 



Defining 



dJ(n) 
do)j 



to be equal to lm(A /n ), A ln equals 
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p-1 



(20) 



To simplify the realization of the above system, 
and to convert it to matrix form, a set of linear filters is 
defined similar to those that were used to derive estimates 
5 for the complex amplitudes of the prominent sinusoidal 
components . 



Z^,z-' =W(e*z) = W(n,z). 



(21) 



Estimates of the frequency deviations Y ln 

(collectively the vector Y n of Figure 1) can be computed 
10 according to: 



(22) 



It is now useful to define the following: 



vfeL-w(n,/Qj 



15 



Thus, 
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>=1 J =-00 

From (23) and (24), we find 



(24) 



(25) 



The vector A n can be computed directly using the 
5 previous estimated complex sinusoidal components {n"}, the 
cross-interference cancelled component estimates X n , and the 
frequency deviation estimates Y n . 

Different methods of applying the frequency 
deviation estimates to determine new frequency estimates may 

10 be applied. In one embodiment the vector A n is used to 
produce new estimates for the frequencies. Referring to 
Figure 1, the AFTs 160 compute the A„ from X n , Y n and the 
previous estimated complex sinusoidal components defined by 
Ip"}. The AFTs then apply the values of A„ to compute new 

15 frequencies. The following equation defines a simple 
gradient method that is preferably used in the present 
embodiment of the invention: 

= & l,n " ^ m &l,n ; 1 = 1. ' * ' . L (26) 

It is to be understood other methods may be used. 
20 In equation (26) // is a very small -sized positive adaptation 
step that is usually a constant value, i.e. it does not 
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vary in time. One may apply the result of (26) into (16) 
and vice versa recursively. However, in the simplest 
implementation, (16) is evaluated and then (26) is 
evaluated, each once at each time instant (step) n. The 
5 overall method is an adaptive method for the estimation and 
tracking of small frequency variations and a least square 
estimation of amplitude of the prominent sinusoidal 
components . 

Referring to Figure 4A and the filters defined by 
10 equation (23) , in some embodiments y ln can be implemented as 
the output of a linear bandpass filter 300 that is tuned to 
the estimated centre frequency a> l , for the I th prominent 
sinusoidal component of the real-valued observed signal. 
Figure 4B shows the substantially equivalent lowpass filter 
15 implementation 310 of the filter shown in Figure 4A. The 
structure of the lowpass filter embodiment is similar to 
that of Figure 3B, but with a different lowpass filter 
impulse response. A third alternative to Figure 3A,4A and 
Figure 3B,4B, is to implement filters in an Intermediate 
20 Frequency (IF) (for both digital and the analogue 
implementation described below) . For analogue 
implementations of the algorithm the IF implementation could 
be advantageous . 

In order to avoid computational error and to 
25 reduce the effect of noise on the generation of the 

estimation of the prominent sinusoidal components QJ 1 = e M " 
the DCOs 140 are used. The DCOs contain recursive complex 
sinusoidal signal generators which operate on previous 
estimated complex sinusoidal components, and the new 
30 frequency estimates a> ln as follows: 
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q; =q;V*'-" (27) 

The output of the DCOs 14 0 is the new set of 
estimated complex sinusoidal components. 

An example DCO implementation is shown in Figure 
5 8. The estimated frequency & ln is converted to a complex 

exponential e J( ° Ln by a non-linear function 800. The previous 
output of the DCO is QJ 1 "' , and this is subject to one unit 
delay 802 so as to be made available at the current 
processing instant. The current complex sinusoidal CI" is 

10 determined by multiplying e J °" by Q"" 1 ' with multiplier 804. 

This structure is repeated for each prominent sinusoidal 
component . 

It is noted that in a conventional FM radio 
receiver, the input signal has only one component and the 

15 amplitude is assumed to be constant. Consequently, a 
structure with behavior similar to a single Frequency 
Deviation Filter (FDF) 120 is used for demodulation of 
frequency. However, the demodulated signal is proportional 
to the deviation of its input frequency and is not exactly 

2 0 equal to the actual frequency. In the approach provided by 
the present invention, this signal is used in a feedback 
loop as in Figure 1 and described by equation (22), to 
adaptively estimate/correct the frequency by minimizing the 
deviation of the observed signal from its estimate. Because 

25 of this, it is expected that the adaptive method for the 
estimation and tracking will provide improvement even in a 
FM radio. 
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Referring now to Figure 2, shown is a flow chart 
that illustrates the adaptive estimation and tracking method 
implemented by the functional blocks of Figure 1. Included 
in the flow chart is a preferred initialization step that 
5 may be used to start the process at its very first 
iteration. 

The adaptive estimation and tracking method has an 
initialization stage 650 and a steady-state operation stage 
660. 

10 The initialization stage 650 begins at step 2-1 

with choosing an incremental step size ja to update the 
instantaneous frequencies of the prominent sinusoidal 
components. The method continues at step 2-2 with the 
selection and initialization of a causal window function 

15 through which to observe the real-valued signal. The window 
function is chosen so that it satisfies the aforementioned 
conditions. An initial estimation of the number of 
prominent sinusoid components comprising the received signal 
is made in step 2-3. This is then followed by step 2-4 in 

20 which the frequencies are initialized. The first three 
steps of the initialization stage 2-1, 2-2 and 2-3 may be 
done in any order; however, step 2-4 can only be done after 
2-3 has been completed. 

The first step of the steady state operation 660 
25 of the method is to observe the real -valued signal, as 
indicated in step 2-5 of Figure 2. In this case it is 
assumed that the iterative process begins with the 
computation of the complex amplitudes. However, it is noted 
that equivalently, the process could begin with the 
30 computation of the frequency estimates. The method 

continues at step 2-6 with the performance of the component 
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extraction step. At step 2-7 the cross-interference 
cancellation step is performed. At step 2-8, the complex 
envelopes are extracted. The output of step 2-8 is a new 
set of amplitude estimates. Next, at step 2-9 frequency 
5 deviation filtering is performed. Step 2-9 occurs on the 
basis of the same set of observed signals x n as was used in 
step 2-6. The method continues at step 2-10 with the 
performance of the adaptive frequency tracking. Finally, 
the iteration finishes at step 2-11 with the updating of the 

10 frequency estimates with the digital control oscillators. 
Depending upon whether or not the entire process is to be 
iterated, as determined by step 2-11, all of steps 2-6 
through 2-11 can be repeated for the same set of observed 
signals. In the simplest implementation, only one iteration 

15 is performed. The method continues back at step 2-5 with 
the observation of the next real-valued signal x n . The 
updating of the complex amplitude estimates is one half- 
iteration, and the updating of the frequencies is another 
half -iteration. 

2 0 Preferably the invention is further enhanced so 

as to be able to simultaneously detect the appearance of new 
prominent sinusoidal components, the presence of already 
identified and previously tracked prominent sinusoidal 
components and the disappearance of prominent sinusoidal 
25 components comprising the real-valued observed signal. 

Referring to Figure 6, the method begins at step 
6-1 with the assumption that at the previous time instant n- 
1, L prominent sinusoidal components have been detected and 
are currently being tracked. The problem is further divided 

3 0 into the following test of hypotheses: 
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1) a previously tracked prominent sinusoidal 
component is still present or has diminished in energy to 
the point where it can no longer be considered a prominent 
component of the real-valued observed signal. In the latter 

5 case, the once-prominent sinusoidal component shall be 
considered to no longer exist within the real -valued 
observed signal and its corresponding parameters will be 
ignored (dropped) and no longer updated; 

2) if a new prominent sinusoidal component has 
10 been detected within the real -valued observed signal, its 

frequency shall be initialized. 

Thus step 6-2 is to compare the estimated energy 
within bands across the periodgram spectrum with a 
threshold. To calculate the periodgram, first over a frame 
15 of signal samples the Fast Fourier Transform of the input 
signal is calculated, and then the periodgram equals the 
squared absolute value of the results. The threshold is 
preferably proportional with <J 2 N (n) . The estimate of noise 
energy, <r 2 N (n) , is obtained by applying the signal 

2 



20 



to the lowpass filter W. , as in (7). This 



allows the identification of all components in the real- 
valued signal which have an energy above the threshold. 

In step 6-3 any newly- identified prominent 
sinusoidal components have their characteristic parameters 
25 initialized within the adaptive estimation and tracking 
method, as previously discussed. 



If the energy of a previously tracked prominent 
sinusoidal component compared with the noise energy, <rl(n) , 
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is very small, then it means that the corresponding 
component is no longer a prominent component of the real- 
valued observed signal and can be ignored by the rest of the 
method provided by the invention as indicated in step 6-4. 

5 An additional step 6-5 can be used in the method 

to improve its performance by estimating the speed of 
frequency variations. The estimated speeds of the 
frequencies can be then used to overcome the crossover 
problem. When the frequencies of two components are not 

10 distant enough their separation may not be possible. 

Consequently, the algorithm may fail to track them properly. 
For each component, therefore, a test may be made to 
determine whether it is far enough apart from other 
frequencies. For each frequency component that is not far 

15 enough apart, instead of using the adaptive algorithm, the 
frequency algorithm can simply assume that the speed of the 
frequency variation is constant during the crossover. 

The five steps 6-1, 6-2, 6-3, 6-4 and 6-5 of 
Figure 6 need to be repeated often enough to keep track of 
2 0 components being added or dropped. The number of signal 
components might change often or not, and as such the rate 
of repetition should be selected depending on the 
application. 

To apply the method to complex observed signals, 
25 the cross interference canceller would be a smaller matrix 
which deals with complex numbers instead, having dimensions 
half of what are required for real valued observed signals. 

In another embodiment, with slight modifications 
the above-described methods/apparatuses can be used to 
30 estimate the chirp of each sinusoidal. Instead of equation 
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(26) above, the following calculations are performed by the 
adaptive frequency tracker 160: 

l,n = (1 " <*)£ /,„-! + " *>/,„-, )l (28) 

where £ /fl is the chirp (the speed of variation of the 
frequency) . In the above frequency estimation procedure co ln 
is estimated as in (26) . Then the variation of the 
estimated frequency in two successive time instant, i.e., 

(^/,n"^/,n-i) is applied to a lowpass filter that provides £ ln as 
an estimate for the chirp parameter of the corresponding 
component. Finally, the frequency for next time iteration 
is predicted by a simple integrator that is the third 
equation in the above procedure. 

In the case that the frequency varies smoothly 
with time the above procedure results in considerable 
15 improvement in frequency tracking at the expense of 

negligible computation. Examples of possible applications 
are wideband FM demodulation, speech frequency tracking, 
chirp estimation in some tracking radars, and several 
biomedical applications . 

Figure 7 is a block diagram of a speech coder 
provided by an embodiment of the invention. In this 
example, an input signal is first converted to digital form 
with analogue-to-digital converter 700 to produce the real- 
valued observed signal x„. This is processed by an adaptive 
estimation and tracking function 702 which operates in 
accordance with one of the previously described embodiments. 
The output of the adaptive estimation and tracking function 
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702 is a set of complex amplitudes 704 and a set of 
estimated frequencies 706. Together, amplitudes 704 and 
frequencies 706 are fed to a signal coding block 708 which 
encodes this information either for storage or transmission. 



discreet signals are being processed and that the entire 
implementation is digital, requiring a D/A converter if the 
original signal is analogue. In another embodiment, the 
apparatus/method is implemented in the analogue domain 
directly eliminating any need for A/D conversion. The block 
diagram of such a system is basically the same as Figure 1, 
except that the component extraction filters and frequency 
deviation filters are continuous time filters having 
continuous impulse responses. Instead of a discrete time 
signal function w k , a continuous time function W t with 
similar properties should be used. Instead of the digital 
controlled oscillators 140, quadrature analogue voltage 
controlled oscillators are employed. 

The adaptive frequency tracker 160 instead of 
employing a delay unit may use a simple integrator as 
following instead of (26) 



Other changes are basically to use dual continuous- time 
elements equivalent to the corresponding digital elements 
that are described above. 

In the embodiments described, it is assumed that 
cross-interference cancellation is performed, and this in 
generally will yield the best performance. In some 



In the embodiments described, it is assumed that 




(29) 
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embodiments of the invention, acceptable results may be 
realized without including cross- interference cancellation. 
For example, in processing signals having very little cross- 
interference, the cross-interference cancellation would 
5 provide only a small improvement. 

In the embodiments described, complex envelope 
extraction is performed at the output to generate new sets 
of complex amplitude estimates. In some embodiments, 
complex envelope extraction can be omitted if a complex 
10 amplitude output is not required. In these embodiments, all 
the steps of complex envelope estimation are performed 
except the last step of extraction, since the output of this 
step is not fed back into other steps /components of the 
method/ system. 

In the embodiments described, it is assumed the 
same weighting function is employed for each of frequency 
tracking and amplitude tracking. More generally, a 
different weighting function may be employed for each of 
these purposes in which case a first weighting function is 
applied for frequency tracking, and a second weighting 
function is applied for amplitude estimation. 

Simulations and Discussion 

Simulations have been conducted for different 
signals. Firstly two scenarios are considered in order to 
study the performance of an example implementation of an 
embodiment of the invention for a signal with known 
components. The first scenario deals with the performance 
in different environments (SNR) , for different signals 
(sinusoids and chirps) and also in the case of a cross-over. 
The second scenario studies the effect of the initialization 
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and investigates the relationship between LIR and the length 
of the window. Then, the example implementation is applied 
to a speech signal as a multi -component stochastic signal, 
to study the performance for speech signals. 

5 Figures 9A and 9B show the results of the first 

scenario where a signal with two components is considered. 
The frequencies of these components are time -dependent with 
a cross -over. This figure shows the performance of the 
algorithm in tracking sinusoids and crossing chirp 

10 components with constant amplitudes for SNR=0 in Figure 9A 
and 2 0dB in Figure 9B . A Hamming window with length 12 9 was 
used. Initial points for frequencies are chosen close 
enough to the true values to ensure initial convergence, and 
the step-size is \i = 10" 4 . The greater the step-size, the 

15 larger the variance of estimated frequencies and the lower 
the bias of the estimation. Figure 9 clearly shows that the 
estimated values converge to the true values. In the chirp 
section, the algorithm tracks the components with a bias in 
the estimated frequencies. This bias/lag is a function of 

20 the window length; the shorter the window, the smaller the 
bias. For the chirp signals a shorter window is preferred, 
while for the sinusoidal signals a longer window works 
better. Around the cross-over moment, the performance 
degrades, since the two components are too close to each 

2 5 other so that they pass through the same filter and cannot 

be separated efficiently. In other words, the CIC is not 
able to cancel out the interference completely as the matrix 
F (Q) becomes ill-conditioned. After the moment of cross- 
over, when the frequencies are far enough apart, the 

3 0 algorithm recovers the frequencies and tracks them again. 
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Figure 10 supports the same results, where the 
estimation errors are plotted for SNR=20dB. The effect of 
lag in the tracking of the chirp components is reflected in 
some bias in frequencies and a jump on the level of the MSE 
5 in both plots. Figure 11 shows the magnitude of the 

estimated amplitudes and the true amplitudes. The estimated 
complex amplitudes are oscillatory with a very low frequency 
resulting from the estimation frequency lag, as they are the 
outputs of demodulators. These oscillations can be 

10 effectively smoothed by allying LPFs to the magnitude of the 
output of the demodulators, assuming a known bandwidth for 
true amplitudes. Figure 12 shows the average of the squared 
frequency estimation error for two different 
implementations. The first one, the same version used in 

15 previous simulations, does the estimation/tracking process 
once in each time step. The second version employs a second 
additional iteration for each time step with a smaller step- 
size of [i 2 = 5 x 10" 5 . The algorithm with a second iteration 
provides a higher accuracy in tracking sinusoidal signals, 

20 as the error variance in frequency estimation is as low as 
10" 7 . Due to the tracking lag for the chirp signals, there 
is a jump in the variance to 10" 5 . Around the moment of 
cross-over, the error increases. After the moment of cross- 
over the error is seen to increase due to the interchange of 

25 the frequencies. Using two iterations per time step, one 
gains shorter convergence period and less lag in chirp 
tracking at the expense of twice the computation. 

Figures 9A, 9B, 10 and 11 show that there is an 
interchange between the components at the moment of cross- 
3 0 over, which the algorithm does not detect. On the other 
hand this interchange is reflected in the estimated 
amplitudes as shown in Figure 11. Hence, using estimated 
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amplitudes the cross-over moment is detectable and 
resolvable as long as the amplitudes are different. 

The second scenario is defined in order to study 
the relationship between the window length and the LIR and 
5 was performed in different environments. Extensive 

simulations show that the LIR does not depend on the SNR 
(LIR measured with a resolution of 0.001/ s where f s is the 
sampling frequency) . For a specific window type, a shorter 
window length causes greater variances in tracking, because 

10 fewer signal samples are. used in the estimation. At the 

same time, a shorter window length provides a wider LIR and 
results in a smaller bias. The tracking bias for the chirp 
component increases for longer windows, because the 
assumption of constant frequencies along the support of 

15 window becomes invalid. 

For a specific window type, since this window 
determines the shape of all involved filters, the LIR is 
inversely proportional to the window length. For instance 
for a Hamming window the LIR is 0.015 and 0.027, when the 

20 window length is 129 and 65, respectively, where the 
frequencies are normalized with respect to / s . If the 
initial frequency error is greater than the LIR, then the 
convergence of the algorithm might take substantial time to 
fall in the LIR. Once the frequency estimation error is 

25 less than the LIR, the algorithm converges with a time 
constant controlled by the algorithm step-size. 

To study the performance of the algorithm for 
speech signals, a Hamming window of length 47 samples is 
used with four components of initial frequencies of 250, 
30 2750, 3500 and 4500 Hz where the sampling frequency was 

11025 Hz. The estimation/tracking process is implemented 
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using two iterations per time step and two step sizes [i ltl and 
|^2,i= M-1,1/3, respectively. Used step-sizes are Hi,i = 0.02 
for the first component and 0.04 for all others. Figures 
13A and 13B show the results obtained for a speech signal. 
5 Figure 13A depicts the tracked frequencies on the spectrum 
background of the speech signal. Clearly, frequencies are 
tracking the dominant energy segments of the spectrum. 
Figure 13B shows the spectrum of the constructed signal and 
supports how successful the algorithm is in decomposing 
10 stochastic signals. 

What has been described is merely illustrative of 
the application of the principles of the invention. Other 
arrangements and methods can be implemented by those skilled 
in the art without departing from the spirit and scope of 
15 the present invention. 



